// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Program  : lin_wacu
// File     : lin_wacu.cpp
// -----------------------------------------------------------------------------
// Revisions:
// Februari 13, 2009  asicherer
// -----------------------------------------------------------------------------
// testbed for the interactive waves model
// -----------------------------------------------------------------------------

// note: as this is testbed-only code: the architecture is not to be taken as an example

#define REMAIN // leave display on after simulation

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <deque>
#include <sys/time.h> // gettimeofday
#include "mermaid.h"
#include "FloatUtils.h"
#include "Array2D.h"
#include "InteractiveWaves.h"
#include "ShipForces.h"
#include "FileReader.h"
#include "defines.h"
#include "PulseControl.h"
#include "Timer.h"
#include "SolverFactory.h"

// colors, shapes and sizes for waypoint/path display
#define PATHPOINT_SIZE    7        // size of path spheres
#define PATHPOINT_COLOR  "orange" 
#define PATHPOINT_SHAPE  "sphere" 
#define STARTPOINT_COLOR "green" 
#define STARTPOINT_SHAPE "capsule"

REAL Fext = 0;

static int steps; ///< number of simulation steps
static int skips; ///< visualizion per 'skips' steps
static REAL currenttime = 0.0; ///< time stamp
static REAL timestep; ///< duration of one step in seconds
static REAL timestddev; ///< standard deviation of time step variation
static REAL pstart; ///< time to start pressure pulse
static int Nx1, Nx2; ///< size in points in x1,x2 direction
static REAL dx1, dx2; ///< grid step size in x1,x2 direction
static REAL grid; ///< minimum of grid step sizes
static REAL g; ///< gravity
static bool running = true; // simulation running?

static InteractiveWaves *pInteractiveWaves = NULL; ///< pointer to waves computing engine
static PulseControl *pPulseControl = NULL; ///< pointer to pressure pulse control


/// reads test file
static bool setup(const char *testfile, const char *solver, Array2D<REAL> **ph0, Array2D<SHORT> **pwet, REAL *ph0_min, REAL *ph0_max)
{
    std::string areafile;
    try
    {
        FileReader fr(testfile);
        do
        {
            std::string section = fr.ReadString();
            if (section == "[Test]")
            {
                steps = fr.ReadInt();
                skips = fr.ReadInt();
                areafile = fr.ReadString();
            }
            else if (section == "[TimeSteps]")
            {
                timestep = fr.ReadFloat();
                timestddev = fr.ReadFloat();
                pstart = fr.ReadFloat();
                pInteractiveWaves = new InteractiveWaves();
                if (!pInteractiveWaves->setup(areafile.c_str(), solver, timestep, pstart))
                {
                    fprintf(stderr, "error setup %s\n", areafile.c_str());
                }
                pInteractiveWaves->getArea(&Nx1, &Nx2, &dx1, &dx2, &g);
                grid = std::min(dx1, dx2);
                pInteractiveWaves->getBathymetry(ph0, ph0_min, ph0_max, pwet);

            }
            else if (section == "[PressurePulse]")
            {
                ASSERT(Nx1 && Nx2, "");
                ASSERT(!pPulseControl, "");
                pPulseControl = new PulseControl(steps, timestep, Nx1, Nx2, dx1, dx2, pstart);
                if (!pPulseControl->setupPulses(fr))
                {
                    fprintf(stderr, "pulses not set up correctly\n");
                }
            }
            else if (section == "[PulsePaths]")
            {
                ASSERT(pPulseControl, "");
                pPulseControl->setupPaths(fr);
            }
        }
        while (!fr.Eof());
        return true;
    }
    catch (...)
    {
        return false;
    }
}


/// test program
int main(int argc, char *argv[])
{
    if (argc < 4) 
    {
        std::string solvers = SolverFactory::getSolvers();
        fprintf(stderr, "lin_wacu testfileName solverName pg|osg|none [capture(0|1)] [realtime (0|1)]\n");
        fprintf(stderr, "solverName may be one of the following: ");
        fprintf(stderr, "%s\n", solvers.c_str());
        exit(-1);
    }
    const char *testfile = argv[1];
    const char *solver = argv[2];
    bool realtime = false;
    bool capture = false;
    capture = argc >= 5 && !strcmp(argv[4], "1");
    realtime = argc == 6 && !strcmp(argv[5], "1");
    if (capture)
    {
        puts("capture");
    }
    if (realtime)
    {
        puts("realtime");
    }

    puts("normal WavesComputer");

#ifdef USE_CUDAINCOMINGWAVES
    puts("with CudaIncomingWaves");
#else
    puts("normal IncomingWaves");
#endif
    
    // setup test and bathymetry data
    Array2D<REAL> *ph0;
    Array2D<SHORT> *pwet;
    REAL h0_min, h0_max;
    if (!setup(testfile, solver, &ph0, &pwet, &h0_min, &h0_max))
    {
        fprintf(stderr, "error setup %s\n", testfile);
    }

    // setup Froude-Krylov module (for now only for first ship)
    mml::Vector3d totalForce;
    mml::Vector3d totalTorque;
    ShipForces *pForces = NULL;
    const ShipGeometry* pShip = NULL;
    if (pPulseControl->getNumberOfPulses())
    {
        pShip = dynamic_cast<const ShipGeometry*>(pPulseControl->getGeometry(0));
        if (pShip)
        {
            pForces = new ShipForces (*pShip);
        }
    }

    // array for elevation data and pressure pulse
    Array2D<REAL> *pPressurePulse;
    Array2D<REAL> *pHeight = 0;
    Array2D<REAL> *pNetHeight = 0;
    Array2D<REAL> *pSpeed = 0;
    int nrp = pPulseControl->getNumberOfPulses();


    // lower pulses automatically if there are only ShipPulses inside grid or MathPulses anywhere
    bool autostart = true;
    for (int i = 0; i < nrp; ++i)
    {
        const ShipGeometry* pShip = dynamic_cast<const ShipGeometry*>(pPulseControl->getGeometry(i));
        if (pShip) // msh ship
        {
            const mml::Vector3d& pos = pShip->getPosition();
            if (pos(0) < 0 || pos(0) > Nx1*dx1 || pos(1) < 0 || pos(1) > Nx2*dx2) // outside grid
            {
                autostart = false;
                break;
            }
        }  
    }
    
    // start timing
    if (realtime)
    {
        getchar();
    }
    int m = 1;
    Timer timer;


            if (autostart)
            {
                pPulseControl->startLowering(0.00000000001);
                autostart = false;
            }

    // Compute the wave evolution
    while (!steps || m <= steps)
    {

        if (running)
        {

            double cgTime = 0.0;
            
            // compute current time
            REAL dt = timestep + timestddev * FloatUtils::randomNumber(true); // new time step
            //    REAL dt = timestep * (FloatUtils::randomNumber(false) < 0.995 ? 1 : 2);
            currenttime = pInteractiveWaves->getTime() + dt;

            // compute pressure pulse
            double moveTime = seconds();
            pPulseControl->movePulses(currenttime);
            moveTime = seconds() - moveTime;
            double pulseTime =seconds();
            pPulseControl->computePulsePattern(currenttime, 0.0, NULL, &pPressurePulse);
            pulseTime = seconds() - pulseTime;

            // compute new wave field
            double waveTime =seconds();
            pInteractiveWaves->compute(dt, *pPressurePulse, &pHeight, &pNetHeight, &pSpeed, &cgTime);
            waveTime = seconds() - waveTime;
//            printf("\r waves:%.1f pulse:%.1f ms ", waveTime * 1000, pulseTime * 1000);            

            //        if (currenttime >= 99.98) exit(0);

            // compute Froude-Krylov force if the first pulse is a panel ship
            if (pForces)
            {
                pForces->computeFroudeKrylov(*pNetHeight, dx1, dx2, &totalForce, &totalTorque);
            }

            // plot the wave surface
            if (m % skips == 0)
            {
                /*REAL */Fext = 0.0;
                for (int i2 = 1; i2 <= Nx2; ++i2)
                {
                    for (int i1 = 1; i1 <= Nx1; ++i1)
                    {
                        Fext = std::max(Fext, (REAL)fabs((*pHeight)[i1][i2]));
                    } // i1
                } // i2
            } // if (m % skips == 0)
            ++m;

        } // running
    } // for all steps or indefinitely


    delete pInteractiveWaves;
    delete pPulseControl;
    delete pForces;

    return 0;
}

